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ABSTRACT 

The published orbits of the planets HD 82943b and HD 82943c correspond to a 
system bound to a catastrophic event in less than 100,000 years. Alternative sets 
of elements and masses, which fit the available observational data and correspond 
to regular motions, are presented in this paper. The planets HD 82943 c,b are 
in a 2/1 mean- motion resonance and are trapped in the neighborhood of a (0,0)- 
apsidal corotation. 

Subject headings: exoplanets, orbit determination, HD 82943, apsidal corotation 

1. Introduction 

This paper is a study of the orbital elements and masses of the planets HD 82943 c,b on 
the basis of the available data. Almost every paper on the dynamics of the actual extra-solar 
planetary systems assumes given sets of elements and masses (which includes the mass of 
the star) as a revealed truth. Sometimes, arbitrary variations of the revealed values are 
introduced in the hope of learning more about the systems stability. Generally, these studies 
are impaired by the scarce information given on the way in which elements and masses were 
obtained. 

The recently published orbits of the planets HD 82943 c,b (Mayor et al. 2004) would 
have a similar fate would they not correspond to a highly unstable system. The published 
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orbits undergo a major catastrophe in about 50,000 years and one of the planets is expelled 
from the system (or collides with the star) in less than twice that time (see fig. 1). In the 
first 10,000 yr of the simulation, the solution lies in the 2:1 resonance zone and one of the 
two critical angles related to this resonance appears in libration; however, the permanence 
of the system in a resonant state is only temporary 

For the sake of comparison with similar simulations, we shall emphasize that the result 
shown in fig. 1 is very sensitive to parameters not fixed by the given data as the epoch 
in which they are valid (which arc arbitrary in a Keplerian orbit determination) or the 
nature of the elements (mean or osculating). It is also sensitive to the actual stellar and 
planetary masses. As it will be discussed later, these masses are only poorly determined. 
The observations allow a function of the stellar and planetary masses and the inclination 
of the orbital plane to be determined, but this function cannot be solved to give separately 
each of these parameters. 

Anyway, in all simulations done, the fate was similar with, at most, a tenfold variation 
in the time of the catastrophic events. Some solutions developed instability in less than 5,000 
yr! However, HD 82943 is a 2.9 billion-yr old star (Mayor et al. 2004) and we should expect 
the planets to be stable for timespans of at least that order. Therefore, no instabihty should 
manifest in a much shorter time scale, be it a backward or a forward simulation. 

Simulations of the primordial disk-planet evolution of pairs of planets in close orbits 
indicate that a capture into resonance followed by a capture into an apsidal corotation is likely 
to occur as a consequence of planetary migration (see Papaloizou 2003). The past evolution 
should be reflected in the orbital characteristics of planetary systems with a period ratio 
not larger than 3 or 4, where the gravitational interaction between the bodies is significant. 
However, the number of systems that could serve as examples of this phenomenon is very 
small. Only two of the possible resonant systems of exoplanets are confirmed: GJ 876 and 
HD 82943. Other candidates, as 47 Uma and 55 Cnc are still affected by doubts concerning 
the very existence of one of the planets (Naef et al. 2004). The small number of such 
systems justifies the need of improving our knowledge on the orbits of the planets HD 82943 
c,b, notwithstanding the several difficulties discussed thereafter, instead of just waiting for 
new elements. A better knowledge of the orbits of these planets is necessary to constrain 
dynamical studies of the evolution of extra-solar planets. 

The results reported in this paper were obtained with a observation series based on 
data published in graphic form. This is a delicate operation and special care was taken 
with points located in the steep branches of the radial velocity curve, since a small error in 
the t-coordinate generates a large O — C {O =observed radial velocity; C = radial velocity 
calculated from the elements) and impairs the results. The radial velocities thus obtained 
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are accurate, but the times of observation are not. It is not possible to get a precision better 
than 1 day from the plot. Anyway, the r.m.s. of the residuals of the series used in this paper, 
with respect to a best-fit solution is the same as that of the original series (6.8 m/s cf. Mayor 
et al. 2004). 



2. The published orbits 

Let us start this paper discussing the orbits published by Mayor et al. (2004). We do 

not need to emphasize again that they correspond to a highly unstable pair of planets and 
are likely untrue. We have to emphasize, however, that the given orbits fit the observations 
very well and that the problem is not that of improving the orbits published by Mayor et al., 
but of finding alternatives to them. As we arc interested in the real planets, we rule out the 
possibility of introducing any blind variation of the elements just to get stable orbits. The 
domain of stable resonant orbits is vast and just pointing out one or more "ad-hoc" stable 
orbits adds nothing to our knowledge of these planets if these orbits do not fit the existing 
observations. 

The orbits given by Mayor et al. come from the fitting of a pair of Keplerian ellipses 
to the observations. In the dual Keplerian model, the mutual attraction of the planets is 
disregarded and the barycentric velocity of the star projected on the line of sight is given by 

Vr = ^ lKfe[cos(/fe + ujk) + ek cosujk] + K- (1) 

k 

where 

nik 27rafe sinik 

The parameters involved in these equations are the system velocity (K), the star mass (M), 
the total system mass (M) and, for each planet, the true anomaly (/), the eccentricity (e), 
the distance from the tangent plane to the celestial sphere (sky plane) to the periastron 
(a;), the inclination of the orbital plane over the sky plane (i), the sidereal period (T), the 
mean distance to the star (a), and the planet mass (m). In the published elements, the 
periastron is measured from the point where the planet orbit crosses the plane tangent to 
the celestial sphere when moving towards the observer (it corresponds to the ascending node 
of the barycentric orbit of the star - near the radial velocity maximum). 

Usually, Kepler's third law is used to eliminate a^. In the case of just one planet, T and 
a are equivalent to their osculating counterparts T, a and we have Air'^a^ — G{M + mjT^. 
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Using this law to separate, in eqn. (2), the known and unknown parts, we obtain. 



sm I 



where M = M + m. The left-hand side: 



(1 - eY'- (3) 



/(M, m,i) = ^ sin^ ^=(--sml) M (4) 

is the so-called mass-function. (See, for instance, Plummer 1916). 

When the system has two planets, and are no longer the osculating orbital elements 
and a correction to Kepler's law is necessary. We have to adopt a different mass-function: 

f{M,mk,ik)^:^sm^ik(M + mk). (5) 

We have, also, to distinguish between the total mass of the system, M, and the sum of 
masses M + and it is no longer possible to cancel the last term in the above equation 
with one power of M as in the case of only one planet. Instead of eqn. (3), we have, 

fiM,m„^,){l + la,) = ^{l-elf'. (6) 

where 0"^ are the classical first-order corrections necessary to transform mean to osculating 
elements (see Appendix). In the case of two planets near the 2/1 resonance, neglecting 
higher-order eccentricity terms, 

(Ti = —0.46— , (72 = 3.0- 



M + rrii M + m2 

We should consider as primary data only those numbers issued from the fit of the dual 
Keplerian model to the observations. They are given in Table 1. Afterwards, they are 
completed with those given in Table 2, which are derived from the primary parameters and 
some additional hypotheses as indicated in the footnote of the table. 

One critical additional datum introduced to obtain Table 2 is the mass of the star. We 
have adopted the mass = 1.15 ± O.OQM© as determined by Mayor et al. (2004) on the 
basis of the Geneva stellar evolutionary models. However, other values exist in the recent 
literature: = 0.93 ± O.OQM© (AUende Prieto & Lambert 1999), I.OSMq (Geneva planet 
search web page, July 31th, 2002) and l.llM© (Chen & Zhao 2002). Several parameters 
given in Table 2 are very sensitive to the adopted stellar mass. Their values are shown for 
two different masses in Table 3 for the sake of giving an idea of the inaccuracies due to 
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Table 1: Parameters determined from the observed radial velocities. (Mayor et al. 2004) 



parameter 


planet c 


planet b 


IK (m/s) 


61.5 ± 1.7 


45.8 ± 1.0 


Sidereal Period (days) 


219.4 ± 0.2 


435.1 ± 1.4 


Periastron time (SJDf ) 


52284 ± 1 


51758 ± 13 


Eccentricity 


0.38 ± 0.01 


0.18 ± 0.04 




124 ± 3 


237 ± 13 


Vr (km/s) 


8.144 ± 0.04 


f SJD (simplified Julian 


date) = JD - 


2,400,000 



Table 2: Parameters derived from those in Table 1. Epoch= SJD 52,000 



parameter 


planet c 


; planet b 


at (10-3) 


-0.7 


4.6 


T (days) (osculating) 


219.2 


437.1 


Mean Longitude at the epoch (deg.) 


18 


76 


a (AU) (osculating) 


0.746 


1.18 


/(M,m,z) (lO-^Mo) 


4.186 


4.113 


msini (Mjup)t 


1.86 


1.85 


msini (lO-^M^)! 


1.54 


1.53 



X assuming i ^ 90° 

t assuming for the star mass M^, = I.I5M0 



the insufficient knowledge of the stellar mass. These differences indicate that the errors in 
the considered parameters are much larger than those that we would obtain by applying 
usual rules to derive them from the errors given in Table 1. It is worth mentioning that the 
osculating semi-major axes and periods, the mass functions and the ratio of two planetary 
masses have only a weak dependence on the chosen star mass. 

For completeness, we should mention that given the small velocity of the system, ~ 8 
km/s, and the short time span of the observations, the apparent and true time scale are not 
significantly different and, thus, we do not need to consider the "planetary" aberration. 
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Table 3: Parameters showing a dependence on the adopted stellar mass 







1.15M0 




1.05Mq 


parameter 


planet c 


planet b 


planet c 


planet b 


a (AU)t 


0.746 


1.18 


0.72 


1.15 


msini (Mj) 


1.86 


1.85 


1.75 


1.74 


msini (10~^M^) 


1.54 


1.53 


1.59 


1.58 



t Osculating 



3. On the 3-body (dynamical) fit 

It is often said that a 3-body (dynamical) fit is significantly better than the fitting 
of the dual Kcplcrian (kincmatical) model to the observations. However, this is not yet 
the case for the planets HD 82943 c,b. This system is much wider than GJ 876 for which 
successful determinations with a 3-body model have been possible (Laughlin & Chambers 
2001). Therefore, the mutual attraction of the planets is much smaller than in GJ 876 and 
the periods are much larger: the observational time span of HD 82943 corresponds to less 
than 4 periods of the outer planet. As a result, the radial velocity curves obtained with the 
dynamical approach or the kinematical one (which correspond to neglect the planet masses 
in the interaction term) are very similar. They are given in Fig. 2 and show that a reasonable 
increase in the observation time span is necessary before getting results different from these 
models. It is worth mentioning that the good observational season for the star HD 82943 is 
the period from November to June and that observations done in the latest season were not 
included in our sample. 

Anyway, we have also run some fits using a coplanar 3-body model. The indetcrmination 
involving M, m, i in the kinematical approach is not completely solved by the dynamical 
approach: the star mass and the orbital inclination continue entangled inside the classical 
gauge Mi, sin i. Since the ratio of the two planetary masses is a robust quantity (not affected 
by the indetermination in the star mass) , it is convenient to substitute one of the planetary 
masses by the mass ratio (Ferraz-Mello et al., in preparation). 

The introduction of the gauge and the remaining planetary mass, as independent un- 
known quantities, in a differential correction procedure has resulted in a correlation factor 
~ 0.99 between the derivatives with respect to the two unknowns. This shows that, with the 
current observations, it is not yet possible an independent determination of the planetary 
masses, exactly as the radial velocity curves of fig. 2 led us to expect. 



- 7- 



Therefore, in the next sections, we will always use the kinematical model, which makes 
the calculations some orders of magnitude faster than the dynamical model and, thus, allows 
a thorough analysis to be more easily done. 

4. Good-fit solutions 

Several authors have insisted on the non-linear nature of the problem under study and 
have used or commented on the possible use of genetic algorithms to determine the minima 
of the function (Laughlin & Chambers 2001; Ford 2004; Gozdziewski & Macicjcwski 2003). 
In the case of the planets HD 82943 c,b, we followed, however, a different path. First, 
we have realized that almost every initial value leads, by a descent along the gradient, 
to the neighborhood of the minimum corresponding to the solution given by Mayor et al. 
Second, the population thus obtained, which should be used to start a genetic algorithm, 
revealed itself good enough to understand the problem. The final population, with 26,000 
solutions, was obtained by means of a random (Monte Carlo) sampling in a wide domain 
of the initial conditions, biased to keep only initial conditions leading to residuals with an 
r.m.s. value not larger than 10.0 m/s. The initial conditions corresponding to r.m.s values 
smaller than 6.9 m/s are shown in fig. 3. This figure presents each solution in the plane 
X2 = e2Cosa;2,y2 = e2sincj2 [^], but it is important to keep in mind that each set of initial 
conditions is complete. Each point in fig. 3 is, thus, the projection on the X2,y2 plane of a 
point in an 11-D space. The variables chosen for this projection, 2:2,1/2, are the variables ill 
constrained by the observations and whose values are critical for the stability of the motion. 
We may compare the distribution of these points with that of the projections of the initial 
conditions on the plane Xi — Ci cosa;i, yi — Ci sina;i shown in the inset on top of fig.3, which 
are much more concentrated. 

The analysis of the results shown in figure 3 involves several complementary aspects. 

First, the [O — C) corresponding to the solutions shown in fig.3 do not differ significantly 
from the (O — C) of the best-fit solution. A straightforward calculation (sec Press et al. 1992, 
chap. 14) allows us to estimate that they belong to the confidence region corresponding to 
the level 38.6 percent - the a/2 level. (The confidence region corresponding to the standard 
1(7 level would include all solutions with r.m.s less than 7.1.) This estimate is founded on 
several hypotheses: a) The distribution of the errors is normal; b) The best-fit solution 
corresponds to the expected value of the quantity — C*)^ (a variable); c) ^{O — CY 
is always the sum of two orthogonal variables. These hypotheses are very stringent and 



The subscripts 1 and 2 refers to the inner and outer planet, respectively. 
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certainly not rigorously fulfilled. 

At this point, we should say that it is culturally difficult to completely abandon the 
heuristics of the least squares. However, we have to keep in mind that our ultimate goal is 
not to look for the minimum of a mathematically defined function, but for initial conditions 
that reproduce equally well the observed values. 

Second, most of the good-fit solutions thus obtained appear grouped in the left lower 
quarter of the plane. We may compare the position of these points and those showing the 
projections on the plane xi — eiCOsa;i,yi = eisina;i (shown in the inset of fig. 3). It is 
well known that anti- aligned periapses and high eccentricities in a 2/1-resonant system are 
ingredients generally leading to instability. Fig.4 shows the fate of solutions whose initial 
conditions are those of solution B (see Table 4), except for the variables X2 = 62 coscj2, y2 = 
62 sin J which are let to vary in the domain shown in figureS. The hachured regions corre- 
spond to orbits bound to collision in less than 260,000 years. The gray strip correspond to 
orbits that survived the simulation time span but that are strongly chaotic and thus bound 
to a catastrophe in a not much larger time span. (The diagnostic of chaos was done with 
the spectral technique discussed by Michtchenko et al. (2002).) A great proportion of the 
solutions shown in Fig. 3, as well as the solution given by Mayor et al., he inside the collision 
region or in the adjacent strip of chaotic solutions. Solutions in the white area behave as 
regular in the considered time span (260,000 years) and we may expect then to survive for 
long times. 



5. Some selected solutions 

We select for analysis two of the solutions shown in fig. 3. The first of them (A) is a 
stable solution taken amongst the solutions with the lowest r.m.s.; it may last for many 
10*^ yrs if the system is edge-on. However, if the mean plane of motion has an inclination 
45 degrees, the planet masses are \/2 larger and the system disrupts in less than 10^ yrs. 
The same occurs if the orbits are nearly orthogonal to the sky plane, but have a mutual 
inclination of 5 degrees. The second selected solution (B) has a slightly higher r.m.s., but 
is very robust and remain regular even if the system has a 30-degree inclination over the 
sky plane, that is if the planet masses are 2 times larger than the minimal masses. These 
solutions are given in Table 4 where elements are grouped by planet to allow us a better 
visualization of the differences in each element. For the same reason some less significant 
digits are kept in the given solutions. 

Notwithstanding some differences in the elements and masses, the O — C residuals 
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Tablc 4: Osculating elements and masses in some chosen solutions compared to those given 
by Mayor et al. Epoch JD 2,452,000 



Solution 


a P e 






m sin i 


r.m.s. 




(AU) (days) 


(deg.) 


(deg.) 


(Mj„p) 


(m/s) 


Planet c (inner) 




A 


0.746 219.44 0.386 


125 


18.2 


1.78 


6.72 


B 


0.746 219.54 0.396 


123 


17.6 


1.71 


6.76 


Mayor et al.f 


0.75 219.4 0.38 


124 


18 


1.85 




Planet b (outer) 


1712 /mi 


A 


1.181 436.90 0.161 


217 


76.2 


1.84 


1.04 


B 


1.180 436.17 0.153 


196 


76.1 


1.82 


1.06 


Mayor et al.f 


1.18 435.1 0.18 


237 


76 


1.84 





Adopted star mass: 1.15 Mq 

t Not osculating. Indicated with M in figs. 4 and 8. 



change only slightly. They are shown in fig. 5 where, for the sake of comparison, the O — C 
taken from Mayor et al. for the same dates are shown. Only a thorough inspection of the 
plots can show differences among them. These figures show that solutions A and B are good 
candidates to be the right orbits and that it is impossible to choose one of them to be labeled 
as the orbits of the planets HD 82943 c,b. The same is true for many other solutions shown 
in fig.3. 

We may proceed and make a study of the selected solutions. The most conspicuous 
results concern the ill determined eccentricity and periastron longitude of the outer planet 
and the critical angles of the 2/1 resonance. The joint variation of the eccentricity of the 
outer planet and the difference of the periastron longitudes of the two planets are shown in 
fig. 6 (the two innermost curves). For comparison we add the curve corresponding to the less 
unstable of the least-squares solutions (r.m.s = 6.70 m/s). The important result we get from 
it is that these solutions are similar and only differ in the amplitude of oscillation of the angle 
Ao; = ujjy — ujc- The virtual center of the loops is a (0,0)-apsidal corotation (aligned periapses) 
with eccentricities Cio ~ 0.13, 620 ~ 0.4, a pair of values consistent with two planets of very 
similar masses (Beauge et al. 2003, 2004). The loops corresponding to solutions A and B do 
not enclose the origin. So, in these cases, the two periapses are oscillating one with respect 
to another, without circulation. 

The curve corresponding to the outer planet in the solution given by Mayor et al. (inner 



-10- 



curve in fig. 7) is very similar to the outer curve in fig.6 (corresponding to a least-squares 
solution). In both cases, the curves enclose the origin. It is clearly seen that the left side 
of the loop crosses the x-axis at the left of the origin. However, the fact that periapses are 
aligned at that moment (the curve passes at the right of the origin) or anti-aligned (the 
curve passes at the left of the origin) has no great signification because the corresponding 
eccentricity is very close to zero (the outer orbit is almost circular). In such case, there are no 
dynamical differences arising from the fact that the periapsis is oscillating or circulating. The 
difference is purely geometric (kinematical) without any topological (dynamical) meaning. 
It is worth mentioning also that when the loop of the outer planet encloses the origin, but 
passing very close to 62 = 0, the loop of the inner planet "leaks" (see fig. 7). It circulates 
in high eccentricity but, again, this is dynamically meaningless because this happens when 
62 ~ and the periapsis of the outer planet is ill defined; the orbit of the outer planet is 
almost a circle and the mathematical location of its periapsis is immaterial. This is certainly 
not the reason for which the least-squares solution is an unstable orbit while our solutions 
are not. 

Table 5: Libration Half-Amphtudes 



Solution 


in ^i(°) 


in ^2(°) 


Period (yrs) 


A 


43 


57 


~620 


B 


72 


89 


-630 



We may draw figures similar to fig. 6 considering the 2/1-resonance critical angles 
9i — 2X2 — Xi — coi, 02 — 2X2 — Xi — 002 [^], instead of Alo. In the 2 solutions, these angles 
librate about 9i — with the half- amplitudes shown in Table 5. Again, when the eccentricity 
of the outer planet becomes very close to zero, the angles 9i and 62 lose dynamical meaning. 
In fact, in all calculations done in this paper the eccentricities and arguments of periapses 
were substituted by the variables Xi,yi introduced in section 4. 

6. Effects due to inclination 

The primary data arising from the kinematical orbit determination do not include the 
plane inclinations (nor the planet masses). However, dynamics simulations depend on the 



^We have written uJi instead of zui as usual in the definition of these variables, since in coplanar models 
the nodal line is arbitrary and the two angles may be considered equal. 



11 



unknown planet masses and the mutual inclination of the planes. The values of these quan- 
tities need to be assumed. For instance, the simulations used to construct fig. 4, were 
done assuming that the orbits are coplanar and i = 90°. In what follows, we will consider 
separately the effects of sin i ^ 1 and the non-coplanarity of the orbits of the two planets. 



If the planet orbits are not edge on, their masses are larger than those assumed to draw 
fig. 4. It is then important to know how the given result change when sini ^ 1. Fig. 8 shows 
two plots constructed in similar way, assuming that the planets are coplanar, but inclined 
with respect to the tangent plane to the celestial sphere. The adopted inclinations in these 
figures are i = 45° and i = 30°, respectively. Therefore, the masses of the planets are larger 
than in the previous case, by factors \/2 and 2, respectively. 

The resonance zones are, now, narrower and conditions for stability are more stringent. 
In the first case {i = 45°), one of the solutions discussed previously (solution A) is close to 
the border of the regularity zone and becomes unstable in relatively short times (~ 10^ yrs)- 
In the second case {i = 30°), this solution is engulfed in the collision region. Solution B 
remains in the regularity zone in both cases. 

A simulation of solutions in the same situation as solution C, with i = 30° (1/ sini = 2) 
showed the same libration pattern and amplitude for longer than 10^ years. 



We have done several simulations with the planets in mutually inclined orbits using 
the RA15 code of Everhart (1985). Mutual inclination contributes to increase the average 
distance between the two planets affecting stability, and we should know how it changes 
the stability of the best-fit solution before looking for alternative fits. The question is of 
how inclined the orbits can be. All cosmogonies discussions are in favor of almost coplanar 
orbits. Anyway, we have explored inclined orbits up to 10° with all possible choices of the 
longitude of nodes. In all experiments, the reference plane was Laplace's invariant plane, 
which is orthogonal to the total angular momentum of the system: 



The order 0(m?) correction of these equations due to the fact that we are using osculating 
Keplerian elements was neglected (For a rigourous formulation, see Michtchenko et al. 



6.1. Coplanar orbits not edge on 



6.2. Mutually inclined orbits 
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2005). The inchnations and nodes are thus chosen such that the component of the angular 
momentum in the reference plane is zero: The two nodes are 180° apart one from the other, 
and 

minia\^l — e\ sinii = 7712/1202 y^l ~ sini2. (8) 

In our simulations, the choice of the individual nodes and periapses was done is such a way 
that, for each planet, the resulting distance uj of the periastron to the sky tangent plane is 
the value given by the orbit determination. Some results posted in Internet, obtained with 
arbitrary node longitudes simply added to the distances cjj, do not correspond to the actual 
planets since the given uji are arguments of periastron only when the reference plane is the 
sky tangent plane. When the arbitrary addition of a longitude of node is done for another 
reference plane, the resulting orbits no longer correspond to the observed ones. 

The results of simulations done with the orbital elements and masses given in tables 1 
and 2 and mutually inclined orbits remained stable for 1-2 million years, but not longer. 
Only one anomalous case exceeded these values and remained stable for the whole duration of 
the simulation (10 Myr). However, it was not possible to reproduce it with a faster computer 
(the same initial conditions and codes leading, then, to a catastrophic event in less than 1.5 
Myr). Besides, in a test of robustness with the same computer, a small change of some 10~^ 
degrees in the position of the reference plane was enough to give rise to instability, again, 
in less than 1.5 Myr. That long lasting anomalous solution, frankly chaotic, seems rather 
due to sticking to an artificial stability island. The exploration of a grid of initial conditions 
with mutual inclinations i < 30° did not show any real island of stability. 

One last consideration concerning the stability of mutually inclined orbits follows from 
the analysis of eqn. (7). First, let us introduce 



£\\ = mifiialy 1 — + m2n2a\^J ^ — (9) 

and write eqn. (7) as 

£||=£ + A£ (10) 

where 

/S.C — 2minia\'J I — e\ sin^ ^ + 2777.2/1202 v/ 1 — e| sin^ 17 > 0. 

It is known that the integral C\\ controls the e-e coupling in a coplanar motion limiting the 
variations of ei and e2, which may remain below some limits (the limit for one eccentricity 
is reached when the other is equal to zero and cannot be exceeded as far as the semi-major 
axes are kept almost constant). If the coplanar motion corresponding to a set of elements 
is unstable because of eccentricities able to reach large values, the existence of a mutual 
inclination acts worsening this condition. Indeed the 3-D integral C controls a similar C\\-i 
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(or e-e-i) coupling (see fig 9). When the inchnation grows, both AC and C\\ grow and improve 
the stabihty of the system (the maximum allowed eccentricities become smaller). But the 
converse is also true. When the inclination decreases, both AC and C\^ decrease allowing 
larger eccentricities to be reached, which may lead the system to become unstable. This 
phenomenon was observed in the simulations with an initial non-zero mutual inclination. 
The system remains stable as far as the mutual inclination remains high, but no bounds 
exist and, sooner or later, it becomes small and the system is bound to a catastrophic event. 

7. Conclusions 

We have shown that several sets of planetary elements and masses exist satisfying the 
observational data of HD 82943 almost as well as the least-squares solutions. The goodness 
of fit of these solutions may be estimated by using some rules, but better than that - 
since we cannot guarantee that the conditions for the use of statistics are fulfilled - is the 
comparison of the O — C plots of these solutions. The proposed solutions are regular and 
correspond to stable planetary systems, at least as far as we may assume that the orbits are 
coplanar and being seen edge-on, at variance with the short lived system that corresponds 
to the published orbits of the planets HD 82943 c,b. If the planets have different inclinations 
and the masses are significantly larger than assumed, solution A becomes unstable and bound 
to a catastrophic close approach of the two planets. 

The analysis of the available observational data shows that this planetary system is 
trapped in a 2/1 mean-motion resonance and the apsidal lines of the two orbits oscillate 
one around another, keeping the system in the neighborhood of a (O-O)-apsidal corotation. 
These orbits suggest a capture into resonance not followed by the complete damping of the 
libration amplitudes. We may compare this system to the GJ 876 planetary system, whose 
hbration amplitudes are close to zero. Either the damping in HD 82943 was not efficient 
enough or it did stop long before a steady state was reached. Since the size of this system is 
larger than that of GJ 876 (the orbits of the outer planet in the two cases are, respectively, 
1.18 and 0.21 AU), we may guess that the larger size affected the relative torques acting on 
the planets and did not allow the system to be damped to a final stationary solution. 

One important question concerns the possibility that new observations allow the real 
orbit to be known. Fig. 10 shows the curves of radial velocities and their differences to the 
curve corresponding to the solution published by Mayor et al. The differences show that 
new observations may help getting best-fit solutions corresponding to non-colliding orbits. 
However, some years will elapse before the observations allow us to distinguish between 
the proposed solutions. Just for completeness of the information given, we mention that the 
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differences were obtained forcing the same period in all orbits to be sure that they correspond 
to real features of the radial velocity curves and not to phase shifts. 

Another point to be stressed here is the separation of the primary elements obtained in 
the orbit determination from those including additional parameters as the mass of the central 
star and the inclination (Section 2). In particular, we mention that the observed periods are 
not osculating and that they should not be used as initial conditions in simulations. The 
recipes to obtain the osculating period and semi-major axes are however simple (as indicated 
in the appendix). 

This investigation benefited from grants from Brazilian Research Council - CNPq - 
Procs. 300575/90-4(SFM), 300953/01-l(CB) and 303196/02-5(TAM), the Sao Paulo State 
Research Foundation - FAPESP - Procs. 98/15431-5(SFM) and 00/07074-0(CB) and from 
Conicet. The authors gratefully acknowledge the support of the Computation Center of the 
University of Sao Paulo (LCCA-USP) for the use of their facihties. 



This is a classical subject in Celestial Mechanics, known since Laplace (at least), which 
has not been considered up to now in papers on extra-solar planetary systems. Besides 
periodic variations, mutual perturbations affect the mean value of osculating elements. Two 
of these variations are particularly important: the mean perturbations in semi-major axis 
and the mean longitude at the epoch. The latter is a secular drift of that angle which directly 
affects the period of the motion. The classical first-order formulas giving these effects are 
the following (Tisserand 1889, chap.20)(see also Ferraz-Mello 1979, sec.5.4): 



where, as already mentioned, a^, Tk arc the osculating semi-major axis and period, is 
the sidereal period and is the mean distance from the star to the planet. Assuming that 
ai < a2 and neglecting higher powers in the eccentricity, we have 



A. Appendix. How mutual perturbations EifFect the periods 
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where a = ^ and b^^if^) is the lowest order Laplace coefficient. For a 0.636 (as in the 
2:1 resonance), we have 65/2 ~ 2.268, dh\i2/da ^ 1.132 and 

(ji = -0.46—— , (72 = 3.0- 



M + mi' ■ M + m2' 
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Fig. 2. — Radial Velocity in kinematical (gray dots) and dynamical (black dots) models 
corresponding to Mayor et al. (2004) data and coinciding at the date 2001.247. The solid line 
shows the difference kinematical— dynamical. The existing observations lay in the interval 
1999.0 - 2003.4. 



-18- 



xl 



-0.25 -0.20 -0.15 -0.10 -0.05 0.00 




-0,25 -0,20 -0.15 -0.10 -0,05 0.00 0.05 0,10 0,15 

x2 



Fig. 3. — Initial conditions generated by means of a biased Monte Carlo technique. The 
main panel shows the projections on the plane X2 = e2Cosuj2,y2 = e2sincj2 of the initial 
conditions whose corresponding solutions have residuals with a r.m.s. less than 6.91 m/s. 
The inset on top of the figure shows the projections of the same solutions on the plane 
xi — ei cosa;i, yi — ei sina;i. 




Fig. 4. — Collision-regularity diagram in the plane X2i 1/2, for the case sin i = 1. The hachured 
region correspond to initial conditions leading to coUision in less than 260,000 yrs. The gray 
region corresponds to solutions surviving that time span but showing strong chaos. In the 
white region, solutions appear as regular. Letters indicate the solutions listed in table tlV. 




Fig. 5. — O — C of the solutions given in Table 4. The crossed squares are the O — C given 
by Mayor et al. The triangles are the O — C of the solutions A and B. 



Fig. 6. — Joint variation of the eccentricity of the outer planet and the difference of the 
periastron longitudes of the two planets in three solutions: the less unstable among the 
discarded least-squares solution (with r.m.s.= 6.70 m/s) and the two selected solutions A 
and B. Polar coordinates: radius vector: 62, polar angle: Au; = Uh — oJc- Prom inside to 
outside: solution B, solution A and the discarded least-squares solution. 



Fig. 7. — Joint variation of the eccentricities and the difference of the periastron longitudes 
of the two planets in the solution given by Mayor et al. (2004). Polar coordinates: radius 
vector: ei (black) and 62 (grey); polar angle: Ao; — cUb — cUc- 
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Fig. 8. — Same as fig. 4 for 1/ sini = -\/2 and l/sini = 2. 
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Fig. 9. — e-e-i coupling during the stable phase of one simulation with initial mutual incli- 
nation 10°. The symbols (alternately black and gray) correspond to the values of the mutual 
inclination. From right to left: mutual inclination ranges 6 — 10°, 10 — 14°, 14—18°, 18 — 22° 
22 - 26°. 
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Fig. 10. — Radial velocity curves corresponding to the solutions discussed in this paper. The 
thick solid lines show the differences of solutions A(gray) and B(black) to the solution given 
by Mayor et al (dotted curve). The wide dots show the dates corresponding to the available 
data. 



